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1. Introduction 

There are at least three areas where standard General Relativity (GR) theory faces 
serious challenges from other competing fundamental theories [I]. Firstly, attempts 
to unite quantum mechanics with GR have so far been unsuccessful and with the 
remarkable success the former has achieved over the last hundred years, many are 
questioning the unique status of GR, suggesting that a fundamental modification of 
the theory could lead to a complete description of gravity on all scales. Secondly, 
in standard particle physics all the forces of nature save gravity have been shown to 
be manifestations of an underlying unified theory of nature and there are hopes that 
modifications of GR could lead to a grand unification of all the known forces into a 
single theory of everything. Thirdly and most important to us here is the missing link 
between the observed universe and the standard theory of cosmology. 

The recent discovery of the accelerated expansion of the Universe has cast a new 
shadow on the simplest cosmology based on Einstein's theory of GR, together with a 
conventional matter description. Observational analyses show that only a tiny 

fraction (~ 4%) of the energy content of the Universe is known to exist in normal 
matter form, whereas ~ 23% of it exists in a little understood dark matter form. 
Dubbed as Dark Energy (DE), the remaining 72%) of the energy content of the 
Universe is believed to be the cause of the inferred accelerated cosmic expansion. 

Many candidates have been put forward as an explanation for DE, but most of 
them fall under one of these three forms: the cosmological constant, exotic scalar 
fields (such as Quintessence) and geometrical dark energy in which the gravitational 
Lagrangian is modified with respect to the usual Einstein- Hilbert one (see [7] , [8] , [9] 
and |10] for an extensive review). 

An important class of modified gravity are the scalar-tensor and f{R) theories 
[TT]-[53]. Both these candidates have their own serious shortcomings |17[ [TBI [T^ . 
and have to pass rigorous theoretical and observational scrutiny before they can 
be accepted as viable theories [24l [25|. In this work we will only concentrate on 
the interpretation of DE as geometrical manifestation of a more fundamental theory, 
focusing on /(i?)-gravity. 

It is well known that the dynamical evolution of small density perturbations, 
seeded in the early Universe, led to the large-scale structure we see today [26]-|36). 
An excellent framework for studying cosmological perturbations is the 1-1-3 covariant 
approach, which has been developed, among other things, to analyze the evolution 
of linear perturbations of Friedmann-Lemaitre-Robertson- Walker (FLRW) models in 
GR [37] - gg. 

In recent years, higher order theories of gravity have attracted a great deal of 
attention. A detailed analysis of the background FLRW models using dynamical 
system techniques has shown that there exist classes of fourth order theories which 
admit a transient decelerated expansion phase during which structure formation can 
take place, followed by a DE-like era which drives the present cosmological acceleration 
(see [17] among many others). However, it has been proved |3H] that when dust 
matter scalar cosmological perturbations are studied in the metric formalism, f{R) 
theories, even mimicking the standard cosmological expansion, provide a different 
matter power spectrum from that predicted by the ACDM model 49 . In [51] 
the evolution of scalar perturbations of FLRW models in fourth order gravity was 
developed for single barotropic fluids using the 1-1-3 covariant approach and the 
solutions of the perturbation equations on large-scales showed that a decelerated phase 
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is not necessarily required to form large scale structures. This divergence from the 
standard GR provides us with a distinguishable signature of fourth order theories, 
which can be tested against observations. 

However, since the Universe consists of a mixture of fluids, a complete treatment of 
perturbations in fourth order theories requires taking this into account. The aim of this 
paper is therefore to present a general framework for studying multi-fluid cosmological 
perturbations with a completely general equation of state in a f{R) theory of gravity, 
using the 1+3 covariant approach. 

This paper has been organized as follows: in section 2 we give the general outline 
of the 1+3 covariant approach. In section 3 we discuss the choice of frame and define 
the key variables used in the description of perturbations in the total fluid and the 
individual fluid components. Equations for these variables are given in section 4. In 
sections 5 and 6, respectively, the scalar and the harmonically decomposed forms of 
our equations are presented. Applications to a radiation-dust cosmological medium 
are given in section 7, with sections 8 and 9 devoted to the analysis of the short and 
long wavelength limits of the perturbation equations. We conclude the paper by giving 
our conclusions in section 10. 

The natural unit conventions {h = c = ks — SttG = 1) are in use. Latin indices of 
tensors run from to 3. The symbols V and ; represent the usual covariant derivative, 
d corresponds to partial differentiation and an over-dot shows differentiation with 
respect to proper time. We use the ( — h ++) spacetime signature in this work. 

For an arbitrary f{R) gravity the generalized Einstein-Hilbert action can be 
written as 

^ ^(/(i?) + 2/:„0- (1) 

A generalization of the Einstein's Field Equations (EFEs) derived by varying this 
action with respect to the metric takes the form 

fGab - t:1; + - Rf) + V,VJ' - gabV^V^f , (2) 

where / ^ /(i?), /' ^ ^ and T™ ^ ^^^^^gf-^- 

Defining the energy momentum tensor of the curvature "fluid" as 



1 



\{f - Rf)9ab + '^b'^af - Qab'^c'^'f 



(3) 



/' 

the field equations ([2]) can be written in a more compact form 

Gab = T^b +Tab= Tab, (4) 

where the effective energy momentum tensor of standard matter is given by 

f^b^^- (5) 

Assuming that the energy-momentum conservation of standard matter T^'^ = 
holds, leads us to conclude that Tab is divergence- free, i.e.. Tab''' = 0, and therefore 
and are not individually conserved [5T| : 



rpm\b ab J_ rpm T^\b rpR;b J rpm r>'-h ( ct\ 

^ab — ji 'JT^^ab^ J ^ ab — 'JT^^ab^ ■ \" ) 
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2. Covariant Decomposition of 4 -order Gravity 

2.1. Preliminaries 

The standard perturbation theory based on the metric formahsm has disadvantages 
when it comes to extracting physical information from the perturbation variables. For 
example, it requires a complete specification of the correspondence between the lumpy, 
perturbed universe and the background spacetime. In other words, this approach is 
gauge- dependent. In this paper we instead use the 1+3-covariant formalism, a fluid 
approach, which, when applied to cosmological perturbations leaves no unphysical 
modes in the evolution of the fluctuations: it requires no prior metric specification 
and is gauge-invariant by construction. 



2.2. Kinematics 

We project onto surfaces orthogonal to the 4-velocity of the fluid flow using the 
projection tensor hab = gab + UaUb and Va ~ h^S/b is the spatially totally projected 
covariant derivative operator orthogonal to m". The covariant convective and spatial 
covariant derivatives on a scalar function X are respectively given by 

X = UaV'X, VaX = hJVbX . (7) 

The geometry of the flow lines is determined by the kinematics of w": 

VbUa = Vfe-Ufc - aaUb , (8) 
VbUa = ^Qhab + (Tab + UJab ■ (9) 

From ([8]) and ([9]) we obtain an important equation relating our key kinematic 
quantities: 

"VbUa ^ -UbUa + l&hab + (Jab + i^ab , (10) 

The RHS of this equation contains the acceleration of the fluid flow iia, expansion Q, 
shear aba and vorticity uJba- 

Another key equation is the propagation equation for the expansion - the 
Raychaudhuri equation (given here for the FLRW background) : 

e + ie2 + i(/i + 3p)-o, (11) 

where fj, and p hold for the total energy density and isotropic pressure respectively. 
This equation together with the equation of state p = p{fi^ s), the energy conservation 
equation 

A + e(/i + p)-o, (12) 

and the Fricdmann equation 

+ ^ - 3a* - , (13) 

form a closed system of equations and completely characterize the kinematics of the 
background cosmological model. 

In this paper angular brackets denote the projection of a tensorial quantity onto 
the tangent 3-space. Thus the relations 

Via) = ha'Vb , (14) 
W^ab) = [hia'h/ - ^h'^hab] W,d (15) 

give the projection of a vector Va and the projected, trace-free part of a tensor Wab 
respectively. 
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3. Matter Description 

3.1. Effective Total Energy- Momentum Tensor 

The thermodynamical description of a relativistic fluid is dictated by the energy 
momentum tensor T^b, the particle flux A^" and the entropy flux S"" of the system. 
Whereas T"'^ and 5"* always satisfy respectively the conservation of 4-momentum and 
the second law of thermodynamics, namely 

T''';6 = , 5";a > , (16) 

particle flux conservation, i.e., N°';a = 0, may be violated. 

The total energy-momentum tensor in a general frame is sourced by /i, p, the 
energy flux g^^t^ , and the anisotropic pressure Tr^^;,^ : 

Tab = M"q"6 + PKb + 2(7(aU6) + T^ab = + T^b ■ (17) 

It defines our thermodynamical quantities: 

= = A„ + /i« , (18) 

ptot^^j^totj^ab^p^^p^^ (19) 

= -n^'h^u^ = C + 9? , (20) 

„tot rptotic j^d ~ni I „R /Ol ^ 

T^ab = ^cd n-ia^) = ^ab + ^^ab i (21) 

where flm = yr, Pm = 9™ = yr, and tt^'^ = yr are the effective 

thermodynamic quantities of matter. 

If we impose the Strong Energy Condition TabV'^V'^ > for all timelike vectors 
then Tab will have a unique unit timelike vector u% {u%Ua = —1). Another 

timelike vector n't can be defined along the flux iV?y, i.e., — , ^ 

" iv ' iv yZ-NbN'' 

For a perfect fluid (or an unperturbed fluid in the background space), uj^ 
and S"" are all parallel |38j and a unique hydrodynamic 4- velocity u° can be defined 
for the fluid flow, in which case 

Tab = piUaUb + phab , N" = nu\ 5" = su" , (22) 

where fi and p are related by the equation of state p — p{fi,s). n — —N'^Ua and 
s = — S'qm" define the particle and entropy densities respectively in the local rest 
frame of an observer attached to u"". 

We can also decompose the EMT with respect to another frame, say but in 
this case we need to introduce a particle drift — h°-bN"' |38[ l40l [52] . 

If the fiuid is imperfect, the fiuid hydrodynamic 4-velocity is no longer unique 
and our EMT will take the more general form given above Eqn. (|17p and the particle 
flux includes a drift term: 

= nu" + f . (23) 

Choosing the relevant frame is a crucial step in the covariant formulation of 
perturbation theories, since u° is the velocity of fundamental observers in the Universe 

§ Fluid flow vector ti° is uniquely defined as the future directed timelike eigenvector of the Ricci 
tensor: = , v^here x''(r) describes the worldline of the fluid in terms of the proper time r. In 
our multi-fluid picture, it corresponds to the normal to the surface of homogeneity. 
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In the particle frame u° = w^, called the Eckart choice, an observer 0„=„„ sees 
no particle drift and hence j"^ — j% — If, on the other hand, we consider the energy 
frame = u^, also known as the Landau choice, an observer O^e measures no energy 
flux {qa = = 0) along the flow line and the EMT takes the forE 

For multi - component matter fluids we have 



ab / , ^ ab J 



(24) 



where 



Tab = M*">fc + PiKb + 1>l + lK + Kb 
Kb = 9ab + U>fc , 



(25) 
(26) 



n,<+j-, (27) 

uj, being the normalized fluid 4-velocity vector for the component, ufu^^ — — 1, 
which we can fix by either choosing the energy frame = u'^^ thereby setting 



0, or the particle frame uf 



u 



Ni 



for which jf 



1^ 

JNi 



for that 



C = TE^ 

component. The velocity of the i^^ fluid component relative to the fundamental 
observer 0„ is defined to be 



(28) 



Vl^ 7^ for tilted, inhomogeneous cosmological media whereas the special case 
where coincides with u°' describes an untilted homogeneous cosmological medium. 
Decomposition of the matter stress-energy tensor with respect to the 4-velocity m° 
gives the following thermodynamical quantities: 

N 



2 = 1 



Pi 



i=l 
N 



la 



(29) 

(30) 

(31) 
(32) 



^Tb = TTdh\A^ = (to first order). 

In a similar way we can decompose the energy momentum tensor of the curvature fluid 
to obtain the corresponding thermodynamical quantities (denoted in what follows by 
a R superscript or subscript). All these quantities, unlike their matter counter-parts, 
vanish in standard GR, with a FLRW geometry: 



R _ rpR a h _ J_ 
M — ab^ 7^ 



^{Rf -.f)-ef"R + f"V^R 



(33) 



P« = Ir^abh'"' = J 



-{f~Rn + f"R + f"'R^ 



I ( Qf'R - /"V^i? - fV^RVaR 



(34) 
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f'RVaR + f'^aR - If'eVaR 
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(35) 



^ab — ^ cd''-{a'''b) — JJ 



f"V^ayb)R + f"yiaRyb)R- CTabR\ ■ (36) 

In the background FLRW universe, V^'^ = and all perfect fluid components have 

Total (resultant) fluid 



Radiation 



CDM 




Figure 1. The Multi-fluid diagram: The different arrows show the unit time-like 
four-velocity vectors at different hypcr-surfaccs SI and S2. The vectors do not 
coincide at the perturbative level. 

the same 4- velocity. By applying the Stewart- Walker Lemma |3S], we can show that 
Vf^ is a first-order gauge-invariant (GI) quantity. If we choose the fluid flow vector 
to coincide with the energy frame (see Fig.l above), then exact FLRW models 
will be characterized by vanishing shear and vorticity of and all spatial gradients 
orthogonal to w° of any scalar quantity [40 : 



0, 



fab — ^ab 

It then follows that, since 

Xa = ^af^ = 0, Ya= V aP, Z, 

in the background, then fi — ^{t), p = p{t) and O 



Vae = o 



(37) 



(38) 



8(t). This necessitates the 
energy momentum tensor having the perfect fluid form, and hence the vanishing of 
the anisotropic pressure iVab and the energy flux qa- 

3.2. Standard Inhomogeneity Variables for the Total Matter 

The key variables characterizing the inhomogeneities of matter are 



D"' = a- 



Ca = a^aR, 
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Aa^VaA, Q^q^-a-^aq, (39) 

where a = a{t) here is the usual FLRW cosmological scale factor. _D™ and Za define 
the comoving fractional density gradient and comoving gradient of the expansion 
respectively and can in principle be measured observationally [40J. The relation 

pe. = + 1 E ^(^'^ - 4)^:^' (40) 

defines the dimensionless variable Sa that quantifies entropy perturbations in the total 
fluid. Wc have defined the shorthand h = /i^ + p^ for the total matter fluid and 
hi = ^i+pi for the component matter fluids, w and denote the effective barotropic 
equation of state and speed of sound of the total matter fluid, respectively, and are 
defined by 

w^^, cl^^, (41) 

whereas for each component mater fluid, these two quantities are given by Wi = pi/fii 
and c% = dpi/d^i. 

3.3. Matter Inhomogeneity Variables for the Components 

The variables characterizing inhomogeneities of matter for the i*''- component fluid 
are defined as 

z,;..^. y;.v.„, 4.«(|:)v...^ (42) 

In near-perfect fluid analyses such as the present one, is often taken to be negligible. 
Thus in subsequent discussions all terms containing this quantity are dropped. 

3.4. Curvature Fluid Variables 

The information about our deviation from standard GR is carried by the following 
dimensionless gradient quantities 

TZa = aWaR , 5Ra = aW aR . (43) 

These variables describe the inhomogeneities in the Ricci scalar. Finally, the velocity 
of the curvature fluid is defined, following [3^, by 

= , (44) 

provided that R 0. Cases with constant scalar configurations or pathologic f{R) 
models with points of inflection in R(t) are excluded from this analysis. 

4. Equations in the Energy Frame 

We now derive the time evolution of the perturbations to linear order in the energy 
frame of the matter, i.e., in the frame where u"^ = n't,. 
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These equations characterize the temporal fluctuations of inhomogeneities in a generic 
perfect cosmological fluid with an equation of state evolving as ti; (1 + w)(w — c^)0. 
They are the following: 



Da + {l+w)Za-wQDa=0 



(45) 



Za -[R■^-'^Q)Za- 



{l + iw)ci - 2(1 + w) fi„i 2cie^ + Scii^iR + Spa) 



1 + w 



2{l + w) f Q{l + w) 



1 Iff" , /V. 



2 /'■ 



&f'{l + w) 



£111 

Rel- 



2K /" 

^7 



1 + w a? 

Ha 



WEa 



f" ~ , Cl^/^D'I" W^^Ea 

'—V^Ua + % , ° + -—^ = 



1+W 



Da 



+ ( 2i?y +e]^a- 



1+W 

3/" 



= 



(46) 
(47) 



1 + w 



-R 



Da 



R Sa 



RZa 



2K 



R'— + R^-' ' ' 



ffff 

7" 



Re- 



Ua-V^Ua^O. (48) 



/" /" 3 /" 

The linearized 3-curvature scalar of the projected metric hab orthogonal to the 
4-velocity vector is [38l 

R = 2(^-^e' + f?j (49) 

reduces to the Ricci scalar in the hypersurfaces orthogonal to when w = 0. The 
covariant, GI gradient Ca gives, to linear order 



Ca 
,2 



f ^ 



£11 £11 

2Qy^a - 2^V27^a 



7^a =0 



(50) 



This variable quantifies the spatial variation in the 3-curvature and is a 
geometrically natural quantity useful in the long wavelength analysis of our 
perturbation equations. The time evolution of this quantity is given by 



Ca=2K- 



3fCa 



a2(2e/' + 3i?/") 



d: 



SwO 



2/^e^ - 6f'fiR 

3(1 + w) 26/' + 37?/" 



6/"V^7^a 

26/' + SRf" 



a\Q /" - 37?/"') ^^R^ff" - 2/" (3/ - 26^/' + 662^^, + 6i?e/" 



3/' 



26/' + 3i?/" 



7^„ 
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66/" . /"I sjj \_^^2 r 36/"7^a mf'DZ 



a2 (26/' + mf") a2 (29/' + 3i?/") _ 

3/"3fia - (e/" - 37^/")7^ J 1 . (5i) 



26/' + 3i?/" /' 

2^2 f 2wa2e ^„ a2 

+ -V < -D H 

3 1(1 + ^'^ /' 

4-2. Component Equations 

These are the equations that describe the evolution dynamics of the individual fluid 
component fluctuations. For the component matter and velocity fluctuations, these 
are given by 

D\ - (m, - cl^)QDl + (1 + w,)Za = -\h,e{cl^iDa + pea) - a{l + w,)Vay''V,' , (52) 

K - [cl - ^) = ^ {h,cl^lDa + h^psa ~ hclfi'Dl) . (53) 

We note that the equations involving the gradients of the inhomogeneities in the 
expansion and curvature variables {Za,TZa,^a,Ca) remain the same as in the total 
fluid equations (|i5|) -([5T |) . This is to be expected since these quantities are global 
intrinsic properties of the spacetime itself rather than of the individual components of 
matter in the fluid. 

4-3. Relative Equations 

Let us now define the variables that relate features of pairs of the different components 
of the fluid, and derive their governing evolution equations. These relative variables 
depend only on the choice of the individual velocities, not on the choice of the overall 
frame. 

Sl^^^-^, in^^in-V^. (54) 

These are the quantities that allow us to distinguish between adiabatic and isothermal 
perturbations [27l l38]. 

The derivation of the evolution equations for the above quantities is 
straightforward and yields 

v:^ - {c% - 1)9^^ - (4 - ct^)ev: = -^(4 - c%h.Dr - -c%si^ , (55) 

SI' + a\7a^'V^' = . (56) 
5. Scalar Equations 

The quantities we have considered so far contain both a scalar and a vector (solenoidal) 
part. Structure formation on cosmological scales is believed to follow spherical 
clustering and therefore we present here the spherically symmetric, scalar density 
perturbations obtained by taking the divergence of the gradient quantities. In so 
doing, we first apply a local decomposition 

aVbXa — Xab — \hahX + S^f^ + , (57) 
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where = X^^^) — ^habX describes shear whereas X[a6] describes the vorticity. 
Vorticity and shear describe the rotation and distortion of the density gradient field, 
respectively. The above decomposition extracts the scalar part of the perturbation 
vectorial gradients. Accordingly, when extracting the scalar contribution the vorticity 
term vanishes 1501. 



5.1. Scalar Variables 

On the basis of the above decomposition scheme, our scalar variables are: 
Z = aV^-Za , C = aV°Ca , 7^ = aV"7^a , 



3f? = aV^SRa , 



The scalar variables describing the total fluid will thus evolve according to 

A™ + (1 + w)z - weA„ ^ , 

2 



z 



■ f" 2 
4-36)^ 



1 1/r , /v» 

2 2 /'2 



/" 



f" - 



(1 + 3w)cl - 2(1 + W) fim 2c2e2 + 3c2(^H + 3pfl) 



2(1 + w) 



/' 



2/'e2 + 3(1 + 3w)^,n + Sf'ifiR + 3pr) 



6/'(l + u;) 



6(1 + w) 
f" 



n 



A, 



1 + w 



-V^A. 



(58) 
(59) 



l + w 

n - 5R + i? 



V^e = , 



w 



l + w 



l + w 



(60) 
(61) 



2R^ + e 1 



+ RZ 

c (a 



(1 - 3c^)a^. 
3/" 



-R 



e + 2^\z 



WUr, 



l + w 

J'" f'ff 



/" l + w 



R 



£ = 



2i?e 



/' 



n 



(62) 



- 2^A™ + 26^3? - 2j^W^TZ = . 
The evolution of the constraint equation is given by 



C = K 



6fC 



a2(2e/' + 3i?/") 



16we 
3(1 + w 



26/' +37?,/" 



(63) 



29/' + SRf" 



2a2(e /" - 37?/"') ^'^RQf.r - 2/" (3/ - 29^/' + 692^;?, + 6i?9/" 



3/' 



29/' + 3Rf" 



n 
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126/" 2/" > I o 

2e/' + 3i?/" /' 



a2(2e/' + 3i?/") a2(2e/' + 3i?/") 



-A,„ + — ^ „ „, 



3(1 + H /' 



3/' 



(64) 



For the scalar variables describing component inhomogeneities and interactions in the 
fluid, the evolution equations are given by 

A;, - Q{w^ - 4)Aj„ + (1 + w^)Z = (c^eA™ + wQe) - a{\ + w^VV, , (65) 

1 + w 



atiti 



,2 ,,i \ i \ 

m J ' 



(66) 



1 



1 



1 



% - ( <, - o ) - (4 - ci,)m = — ^(ci - 4-)m.a:„ - -4.5,, , (67) 



5,, + aW^Vn = 



(68) 



5.2. Second- order Equations 

The above first order equations (|63l) - ([55)) can be reduced to a set of linearly 
independent second order equations. This has the advantage of simplifying the 
equations and making comparisons to GR more transparent [51) : 



(c2 + ^ - 2w)e - 



A - r^V^A 



f" 12K\ 



/' J 



1 + w 



f" 



l + {f-2tim + 2Ref")jj^~2Re[^] -2RQ'— 



/' 



/' 



n 



(1 + - (27^ + f - 7 - 5 ) + ' 



(69) 



n +[2R'— + Q\n+[R^— + ie^ ' "^'^ ' ■' 



) J 

f" 

c? - 1 . - 

-RAr, 



/" ■ " /" ■ r ^3/" 3 . 



1 + w 
R 



(3c2-l)/i„. , ^ + ^lj^Q_, 



3/" 



1 + w 

W^jri , 2w - cl ■ 



A, 



/" 1 + w 



w., QA, 



RQ 



1 + w 
I A,„ - V27^ 
2w 



1 + w 



-Re 



R + R^-' 



f" 



1 + w 



rill 



1 + ?« 

e = 



1+Wi 

1 + w 



f" ~ 

Am-{l+W,)yV^n 



(70) 
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/ 



(i + -)^-(2^--^-2ei?V)c^ + c?e 



ifr_ 

2 /'2 



/'• 



£111 

ReL. 



7^ 



(w - 



-c,,w)Q ~w[2—- — -m — 



£ = 



(71) 



The second order equations (|69p - (l7ip governing the propagation of the entropy 
perturbations for a general e (or Sij) are in general very comphcated and consequently 
we will present them for specific (radiation-dust) applications in section [T] 



6. Harmonic Analysis 

The above evolution equations can be thought of as a coupled system of harmonic 
oscillators of the form 

X + AX + BX = C{Y, Y) , (72) 

where the second term from the left represents the friction (damping) term, the third 
one, the restoring force term while C represents the source forcing term. A key 
assumption in the analysis of the equation here is that we can apply the separation of 
variables technique 

X{x, t) = X{x)X{t) , y(x, t) = Y{x)Y{t) , (73) 

and write 

A = ^A'=(t)Qfe(f), Y ^Y.^\t)Qk{S) , (74) 

k k 

where Qk{x) are the eigenfunctions of the covariantly defined Laplace-Beltrami 
operator on an almost FLRW space-time: 



(75) 



Here k — is the order of the harmonic and Qk{x) = (Q is covariantly constant). 
In this way the evolution equations and the constraint equation can be converted into 
ordinary differential equations for each mode. After harmonic decomposition the first 
order total and component fiuid equations (|63p - (p5|) can be rewritten in the following 
form: 



A^, + (l + ?«)Z*-u;eA,^ 







(76) 



R 



f" 



e z" 



/' 3 
(1 + 3?i7)c2 - 2(1 + w) fl„ 



2(1 + w) 



/' 



2e^9^ + 3c^(Mi? + 3pfl) 
6{l + w) 



1 + w a'^ 



2fe^ + 3(1 + 3w)^irn +3f' inn + 3pr) 1 P 
6/'(l + u)) 1 + 



we 



f" 
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1 fcV" 1//" , />™ U(,(f"\\up>f"' 



VJ' = 0, 



1 + w 



-A" 



= , 



1 + 7« 

(1 - 3cl)fi 



(77) 
(78) 
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1 + ?« 
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A'' 



1+W 
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/ 1.2 fill f{iv) 



RZ" + I ^+ R— + R 
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(79) 
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a^{2Qf' + 3Rf") 
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129/" 2/" 
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16w9 4/'92 - 12/Vi? 



3(1 + w) 29/' + 3i?/" 



12/" 



29/' + 3i?/" 



27^^a2(9/" - 3i?/'") 12^e/'/'" - 2/" (s/ - 2/'(92 - 3^7^) + 6i?9/" 



fc2 

^2 



3/' 
2>{l + w)' 



29/' + 3i?/" 



2a2/".jjjfe 2a2(9/"-37?/") 



3/' 



(81) 



1 + w, 



Af - (u;, - 4)9A,^ + (1 + Wi)Z'' = {clAt + we'') 9 + (1 + Wi) — V, , (82) 

1 + w ^ 



wfe „ ^v-/? = 



(83) 
(84) 
(85) 



The form and use of Eqn. (|8ip will be more transparent when we discuss the long 
wavelength limits of our perturbations for radiation and dust backgrounds. The 
harmonically decomposed set of second-order equations (l5^ - d7T|) will become 



A* 



{cl + \- 2w)e - R^ 



2 1 r 2 A 1 \ 

w + 5c„ — 4w — 1 ' 



/' 



._(3y;_5c2)^ + ,c- 



2 , , • /" 12X\ 

2 ^)i2R~AReL--^j+c^- 
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+ (c^ - 4) e 



A''^ + (l + wOe^7t'= 
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As can be seen, this second order set of equations is not closed. For a two-component 
fluid the entropy and velocity perturbations equations are given by 

fc2 



fr2 

Qfe - —V ■ 



(89) 



V/' 



c2. - r2. 

SI SJ 



1 



' a(l + w) V3 

1 
3 



w - GA, 



A I ^2 9 3c2 



a 



c2. - c2. 
a[l + w)" 



A. 



(90) 

Since Eqns. ([M)) and (jHOl) are not linearly independent equations, we can choose 
either one of them to close our system of second order equations (|86ll88p . 
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7. Perturbations in a Radiation - Dust Universe 

7.1. Background Setup 

Now that we have derived the equations for perturbations of a general muhi-fluid 
system, we consider an apphcation of the equations for a cosmological medium 
containing a non-interacting radiation-dust mixture and described by a flat [K — 0) 
FLRW spacetime. Since our component fluids satisfy the conservation equations 
separately, we write 

Ad + e^id = , (91) 

iir + ^6/^^ = , (92) 

where d and r subindices hold for dust and radiation respectively. 

The general equation of state w for such a radiation-dust mixture is given by 

Pm Pd+Pr 1 Mr 
W = = = (93) 

Mm IJ'd + ^J-r 3 Md + Mr 

and the speed of sound in the mixture is 

cl^tnL^ fMr ^ j-g^^ 

M„ 3(3Md + 4/Ltr.) 
Wherever necessary, we will use the shorthand 

\ {hrcla + hdcl^) . (95) 

In general, since we do not have an explicit expression of the Hubble parameter 
H and the curvature scalar R as functions of the scale factor a in generic f{R) gravity 
theories, an exact multi- fluid background solution is not available and numerical 
solutions need to be obtained. This important issue will be investigated in a future 
work. 

In this paper, we will confine our discussion to i?" models [HI [5T] and 
look for solutions in the short wavelength and long wavelength approximations for 
perturbations deep in the radiation and dust dominated epochs. During these epochs, 
since one fluid is negligible with respect to the other, we can use the exact single fluid 
background transient solution for i?" models given by a = aeg(t/te(}) where Oeq 
is the scale factor at the time of radiation-dust equality teq and will henceforth be 
normalized to unity. 

In i?" models the expressions for the expansion, the Ricci scalar, the curvature 
fluid energy density, the curvature fluid pressure and the effective matter energy 
density are given by 
2n 

(1 + w)t 

_ 4n[4n-3(l + H] .97. 
2(n~l)[2n(3^;-f5)-3(l + H] .Q„^ 
2(n-l) fn(6u;2+8u;-2)-3w(l+w)] , , 

PR = on , ' 99) 

3(1 -I- wyt'^ 
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- 3n(l + w) 

(1 + W)2t2 



4n3 - 2n{n - 1) [2n(3ui + 5) - 3(1 + w)] 
3(l + w)2t2 



.(100) 



7.^. To to/ Fluid Equations 

Upon expanding Eqn. (|40p for a mixture of dust and radiation, we obtain 



3(3^d + 4^r) 
and hence 



e = — - 



-S, 



dr 



Sdr , (101) 

(102) 

3(3/id + AflrY + 4A<r 

Using these relations and applying the general total fluid second order equations 
to the radiation-dust mixture yields 



3^d + 4^r 

We can thus readily derive the evolution equation for e as follows 
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c^RS — , 



s. 
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dr 



1 



^dr 



Ink 
2 ^z^dr 



where A„i and Sdr are given by A^ = ^""^l^^l^^^'' , -^dr = A^ - fA^ . 



(105) 
(106) 
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The perturbations of the density gradients of the radiation component of the fluid are 
described by the propagation equation 



A k 



1 _ 4 / ^wfid 
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Similarly the propagation equation of the dust density gradient is given by 



2 jid , 4w 

3 l + «;4/id + 4Mr 



1 



+ 



l + w 

Hd 
l + w 



V3Aid + 4/Ur lld + IJ-T 



e + 



) e. 



fj-d 



R 



f" 



(l + «;)(Md + Mr) /' 



( w 



Aw 



■ f"' 
-R— 

fJ-d + Mr /' . 



\Md + Mr 3Md + 4/Ur- 3(/id + /ir 



R& 



/" 

— H 



4Mr 



Awnr + {Aw + l)fid 4(c2 - w) (c^ - w)c2 



Md + Mr 

1 + 7« - 2c2 
Md + Mr 



3Md + 4Mr 



Md + Mr 3(Md + Mr)^ 



3(3Md + 4Mr)2 



8w 



Slid + 4Mr y /' 



fJ.d + /Jt 



Aw 



I 

3Md + 4Mr IJ'd + IJ'r J f 



a5 



Covariant gauge-invariant perturbations in multifluid f{R) gravity 



1 



1 + w 



(w - cDur 



lid + Mr + 3(/id + llrY 



■ f" 



3(Md + Mr)^ 



Slid + 



I 



AfldlJ-r (1 - 3w)Mr - Swfld 
{Slid + 4:IJt)^ 3(/id + Mr) 

' (l + w - 2c^) Qwnd \ 



IJ-d + IJ-T 



f 



1 , k^f" Iff" , /"(Mr + Md) 



+Re 



r 
f 



2 a2 /' 2 /'2 

7^'= + e^Tt*^ = . 



+ 



f,2 



RO ( 47 



/' 



Covariant gauge-invariant perturbations in multifluid f{R) gravity 



20 



In terms of the component perturbation variables of section [3] we can rewrite the 
propagation equation for the curvature fluid gradient as 



fc2 .. /"' . 2 /(») . /"' 1 /' R 
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8. Short Wavelength Solutions 



In this section, we will study the evolution of the short wavelength modes, i.e., large 
values of the wave number A:, by using the equations presented in section 6 valid for a 
radiation and dust mixture. The general results will then be considered for i?" models 
and a proposal for a quasi-static approximation for the matter perturbations will be 
introduced for both radiation and dust dominated epochs. In that approximation, 
widely used in the literature [42], all the time derivative terms for the gravitational 
potentials are discarded, and only those including density perturbations are kept. 
The decoupling process for the involved equations is therefore highly simplified. 
Nonetheless, when this approximation was used to study fourth order gravity theories 
in the metric formalism, it was proved as too aggressive and a more detailed analysis 
is required [48ll43] . 



8.1. Perturbations in the Radiation-dominated Epoch 

Let us now look at the case where the characteristic size of the fluid inhomogeneities 
is much less than the Jeans length for the radiation fluid but is still larger than the 
mean free path of the photon, i.e., A <C Xh <C Aj. Similar investigation has been 
made by [41] for the case of GR. Note that the scale dependence of the perturbations 
equations is not trivial as can be seen in |44j . 

Here we assume that we can neglect the interaction between the component fluids 
and that the radiation energy density can be taken as almost homogeneous, meaning 
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w 0. 

This amounts to studying dust and curvature fluctuations on a homogeneous 
radiation background, whereby radiation affects the growth of the dust fluctuations 
by speeding up the cosmic expansion [3^. We consider the flat {K = 0) background, 
and hence the equations for such a background are given by 
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In these Hmits the above set of equations (|110mi5|) can be rewritten as 
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(121) 
(122) 

(123) 



From Eqns. 
equations: 



/" /" /" 

(|118p - (|123p we obtain the following two second order differential 
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7^''■ = 



(125) 



/" /" /" 3/" 

It can be shown that H and Rf'/f are of the same behaviour for i?" models, 
whereas fid fJ-n implying that curvature and radiation fluids effectively dominate the 
fluctuation dynamics. In effect, terms like /i^A^ merely sub-dominate in the curvature- 
radiation-dust mixture. Hence we can safely approximate the above equations by 
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(126) 



= 



(127) 



These two equations tell us that, deep in the radiation-dominated era, the curvature 
fluctuations are decoupled from the matter in the second order equations. 

GR is a specific example of the generalized i?" models where n = 1. In this limit, 
Eqns. ([T261) and (fT26| reduce to 



1 



+ 2i/A^ + _7^^ = , 



(128) 



7^'= = , (129) 

thus yielding the standard GR equation for the density contrast in a radiation 
background 

A^ + Ia^ = , (130) 
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whose general solution is given by 

A^(t) =Ci+C2lnt. (131) 

with Ci^2 arbitrary constants. For n 7^ 1, with it; = i in the radiation-dominated 
epoch, Eqns. (|126p and (|126p take the following forms: 
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(132) 
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7^'= = , (133) 



.2 

where we have used the fact that ^ - 
the time of radiation-matter equality. 



t " with normalized time t„ 



1 at 



Quasi-static Analysis In general the system of equations (jl32p - (jl33p yield Bessel 
hypergeometric type analytic solutions. However, since we are dealing with small 
scales we can take a quasi-static approximation where the time variations in H!' are 
treated as negligible, i.e., TV' ~ and TZ^ ~ 0. In this scheme the overall dynamics of 
the density perturbations leads to the simplified, k- independent, equation 



(134) 



This equation admits the general solution 

A^(t) =Ci+ C2ii-" . (135) 

On small scales, radiation suppresses the growth of fluctuations as they enter the 
horizon before radiation-dust equality, and dust (baryon) self-gravitation is not yet 
strong enough to offset the cosmic expansion. This is because the expansion scale 
factor grows faster than the perturbation amplitudes do. The phenomenon is known 
in the literature as the Meszdros effect [53] . 

It is clear from the above analysis that the Meszaros effect puts a constraint on 
the value of n in i?" gravity. To do so, all we need do is determine the allowed values 
of n for which the perturbation amplitudes grow slower than the expansion in the 
radiation dominated era, i.e.. 



d 
dt 



ait) 



dt 
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3n 

< ^ 1 < 

2 



(136) 



This means that only values of n > | give a growth rate compatible with the Meszaros 
effect. 

In figure 2, we plot the normalized dust density contrast S{t) 
radiation-dominated epoch. 



A™(t)/Ae5 in the 



8.2. Perturbations in the Dust- dominated Epoch 

During this epoch of the universe the dust energy density is dominating in the two-fiuid 
dynamics and all order-of-magnitude approximations go in line with the assumption 
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Figure 2. Dust growth factor in the radiation dominated epoch for R" 
models: The plots show the growth factor obtained by solving numerically the 
full system of equations I I132I I and II133I I for scale A^f = lOOA and the quasi-static 
solution l|134| l for n = 1.5, 1.4, 1.3, 1.2, 1.1 from bottom to top. The top-most plot 
corresponds to GR (n = 1). It can be seen that quasi-static results are quite 
close to those of the full system for the stated values of n, only slightly (but 
invisibly) lower in the plots. For values of A^f > lOOA the growth factor appears 
to be insensitive to scale showing the convenience of introducing the quasi-static 
approximation. 



that fjid » fJ-r- The equations (|110|) - (|115p will, upon imposing the short-wavelength 
assumptions (|116lll7p . become 
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As can be observed, these two equations differ from their counterparts in the radiation- 
dominated epoch in that the curvature perturbations are not decoupled from that of 
matter in the system of equations. 

The limiting GR perturbation equations for (|143p and (|144p in this epoch are 
given by 
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This equation admits the well known solution 

A^(t) ^Cit-^ + C2t^ . 
For i?" models, equations (|143ll44p take the form 
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Quasi-static Analysis In the quasi-static limit with (-f >> 1 we get a single 
second order fc-scale independent equation 



4(8n2 - 13n + 3) 



A;5 = 0, (151) 



9^2 

the solution of which is given by 

AS(0 +C2i"- , (152) 

where a± = + i ± V-ii2«^+i84n ^^ rj.^^ coefficients Ci,2 can be determined by 
imposing initial conditions. 

At t = teq = 1 we have 

Afrf)e,^Af^)(ie,)=Ci+C2 , (153) 

and differentiating (|152l) gives 

AS(t) = a+Cir+-i + a_C72i"-~' , (154) 
which, at equality, will give 

Af,),^ ^ Afrf)(ie5) - + a_C2 . (155) 

Solving (|153p and (I155P simultaneously we obtain 

C, 2 = . (156) 

Fig. 3 shows the evolution of the density perturbations 5{t) = A^^^{t) / A^^^^^ in time 
(t from 1 onwards, where t = t^q — I is the normalized time at equality) for the above 
linearly independent solutions, Ci_2 having been obtained by setting A^^^^^ = 10~^ 
and A^,)^^ = 10"^ . 

9. Long Wavelength Solutions 

For specific intervals of n, a set of initial conditions give rise to cosmic histories which 
include a transient decelerated phase which evolves towards an accelerated phase. 
Structure formation takes place during the transient regime [SJ . 

In this section we analyze the evolution of scalar perturbations during this 
phase, in the long wavelength limit. In this limit the wavenumber k is so small that 
A = 2i£ ^ Xh, i.e., ^^fji ^ 1- AH Laplacian terms can therefore be neglected and 
spatially fiat {K ~ 0) backgrounds guarantee the conservation of C, i.e., C'^ = 0. In 
this paper we are considering only adiabatic perturbations, i.e. Sij — and hence, for 
a radiation-dust mixture, the equation for the evolution of entropy perturbations 

Sdr + aV^Vdr = . (157) 

implies that 

Vdr^O. (158) 
And from this and the equation 

Vdr - (cl - \) eVdr = -\{cld - cDtlA^ - -clSdr ■ (159) 

\ 3/ ari a 

follows 

{eld - cDmA™ = . (160) 
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Figure 3. Dust growth factor in the dust dominated epoch for R" models: 
The plots show the growth factor obtained by solving numerically the full system 
equations I I149I1 and I I150I I for scale A^f = lOOA and the quasi-static analytic 
solution ijlSll l . It can be seen that quasi-static results are indistinguishable from 
the full results for n = 1.5,1.4,1.3,1.2,1.1 from bottom to top, with the full 
system solutions slightly higher than those of the quasi-static approximation. It 
can also be seen that for higher values of n the growth factor increases more slowly 
till a critical value of n somewhere between 1.4 &: 1.5 where the growth factor 
becomes a decreasing function of time. Note the n = 1 case (GR) is presented on 
the topmost plot. 



We therefore have the following system of equations: 
A™ (1 + w)Z - wOA™ = , 



Z-{R^j-'^@\Z- 



{2cl ~W-\) Hr. 



f" 



l + w 



1 Iff" /V 

2 2 /'2 /'2 



i?A,„ = , 



(l+w) f {l + w)\2 f 

2 



/' 



-Re 



Re 



r 
f 



7^ = 



3/" 



iRe + R)Lj, + R'L_ + ^-j 



n = o 



(161) 



A, 



(162) 
(163) 
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Co ^Rf" 



2^e^-^(-^-2^ + 2i?e^ 



7^ 



(164) 



Co being the conserved value for the quantity C. 

In terms of the background i?" solutions and making use of the conservation of 
C the above equations can be rewritten as 



1 + w — 2n 
1 + w 



6(n — l)r 



n + 3(n — l)^ — 3 

3(l + ii;)2 



t 



4al [n + 3(n - l)w - 3] [An - 3(1 + w)] 



1 4n 



9(n- l)(l + u>)3i 



3+2 



+ 



7e = 5R + 

5ft = 
+ 



4 [n + 3(n - 1)m; - 3] [4n - 3(1 + w)] 

3(n - 1)(1 + wf [n{6w + 8) - 15(1 + w)] 
4 [n + 3(n - l)w - 3] [4n - 3(1 + w)] 

8ncl [4n - 3(1 + w)] 

3(l + w;)3 ' 

(n - 4) + 2(;) - 2)w 3;)(;(-l) 



tn . 



(165) 
(166) 



(l + w)t 

2n(4n - 3w; - 3) 



3? 



- 2 



(1 + w)[n + 3(n - l)w - 3] 
9n(n-2)(n- 1) 



n + 3w{n — 1) — 3_ 

4n o 

Cot~^^^+^ 

3^2 (9n - 26) + 57n 



n + 3(n — l)w — 3 



+ 2n^ - 7n - 



8n2(n-2) 



9(l + w)(n-l) 9(l + w;)2(n- 1) 



7e 



A„ [4n - 3(1 + w)] [An + 3(n - l)w - 3] 



(^n - 1)(1 + w)^ [n + 3(n - l)w - 3] 
[(9w;(l + w)+8)n^ - {27w^ + 24w + 13) n + 3(1 + w){l + 6w)] . 

9.1. Perturbations in the Radiation- dominated Epoch 



(167) 



The second order set of equations governing the dynamics of density perturbations in 
this epoch is given by 



A, 



, n(9n - 14) + 4 . n [n(n(19n - 54) + 58) - 32 ] + 8 ^ ^ 

2(n - 2)2f2 



+ 



2(n- 2)i 

2 [n(3n - 4) + 2] 
3(n-2)2 



-A" 



_ n(15.-22) + 14^, ^ lfci).-"Co = 0, 



3(n- 2) 



3(n-2)2 



(168) 



_ n(lln- 32) + 32 ^fc ^ 3 [n(5n - 9) + 8] ^^ _ 3n - 3) + 2] 



2(n-2)i 



2*2 



2(n - 2)^3 
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Making use of the eonservation of C, we can eliminate TZ'^ and TZ'^ quantities in favour 
of (and its derivatives) and Co- This way we can get a decoupled third order 
k-scale independent equation for A^: 

n - 5 £ 24n - 19n^ + 8 . - 2) [5n^ - 8n + 2] (12 - 7n)Co 



This equation admits the general solution 

A';{t) = Cit^-^ + C2t^+ + Cat^- + C\t^-" . (171) 

where Ci,2,3 are arbitrary integration constants to be evaluated from initial conditions 
with 

_ 2(24-14n)Co 

and 

^^^^^_^n^ym^^^±m. (173) 

Provided that the initial values of A^, A^, Ajr and Co axe known at tgq = 1, the 
integration constants can be determined since 

A(r)e« = Cl+C2+C3 + C4, 



K)e, = [ '-Y^ j + C2/3+ + C3/3- + (2 - n)Ci , 
(n-2)(n-4)" 



n - 2 



(»')e9 



Ci + C2/3+(/3+-l) 



+ C3/3-(/3- - 1) + (2-n)(l -n)C4 . (174) 
We do not present Ci,2,3 explicitly for the sake of simplicity. 

9.2. Perturbations in the Dust- dominated Epoch 

Proceeding in a similar fashion for the dust dominated, long wavelength, regime gives 
the second order evolution equations given by 
■^^ , n(8n-13) + 3 ^, , K8n - 13) + 3] Kl6n - 15) + 9] ^ ^ 
+ {n-3)t + 3(n-3)2f2 

3(n - 1) [n(16n - 15) + 9] ^j, _ ^^ [(n(16n(8n - 31) + 711) - 540] + 189 ^ 
^ 4(n - 3)2(4n - 3) * 4(n - 3)2(4n - 3) 



n (27 + 54n - 5671^) - 27 
4(n - 3)2(4n - 3) 



t-S-Co = , (175) 



_ 4(n-l)[n(2n-5)+6] -fc 4 [n (n(2n(16n - 65) + 213) - 198) + 81] , 
[n(n - 4) + 3) t 9 [n(n - 4) + 3] 

16n(3-4n)2[n(8n-13)+3] ^fe 2n [n(4n - 7) + 3] ^_(„^2)^ . 

27[(n-4)n + 3]f4 " n(n - 4) + 3 * ^° " ^ ' ^^^^^ 



(170) 
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which reduce to a single third order evolution equation for the density perturbations 
given as 

^Afe + l^A^ - 2[^(4n(8'^^19) + 33) + 9] d , _ 2(4n - 3) K8n - 13) + 3] ^ 
dt^ t df^ 9(71 - l)i2 dt 9(n - l)i3 

(I2n2-31n+18) , , 

- -^TT 7^-—r^Ca = , 177 

6(n-l)t"+i ^ ^ 

which is a third order decoupled k-scale independent equation. The general solution 
of (|177p is given by 

A^(t) = Cit-^ + C2t^+ + Cst^- + Cit^-"^ , (178) 

where Ci_2.3 are arbitrary integration constants to be evaluated from initial conditions 
and 

^ 9(l2n2-31n+18)Co 

8(48n4 - 184n3 + 159n2 + 63n - 81) ^ ' 

together with 



1 1 /256n3 - 608^2 + 417n- 81 

^±"-2^6V ^^1 ■ (^'^^ 

As in the radiation epoch, the integration constants Ci^2,3 can be determined 
from the initial values of A^, A^, A|; and Cq known at t^-q = 1 as follows: 

^\d)eq = Ci + C2 + C3 + C4, 



^Ue, = -Cl + + C37- + (%^) C4, 

Af.w = 2Ci + C27+(7+-l) 



+ C-»-(7--l)+ "'--'"f-''"' ft. (181) 

Once again, for the sake of simplicity, we do not present them here explicitly. 

It turns out that in the adiabatic limit, the long wavelength solutions of the 
growth factor both in the radiation and dust epochs are exactly the same as those 
found in |51| . 

10. Conclusions 

In this work we have for the first time presented a detailed analysis of the (1 + 3)- 
covariant and gauge-invariant theory of cosmological perturbations in situations where 
the universe is described by a multi-component fluid, with a general equation of 
state parameter for an arbitrary /(i?) theory of gravity. The linearized evolution 
equations of the density and curvature perturbations of such a universe have been 
derived for both the fluid components and the total matter, relative to the energy 
frame. We then have taken the background transient solutions of i?" gravity for a 
two-fluid system dominated respectively by radiation and CDM (dust) and obtained 
solutions in both the short and long wavelength approximations. These solutions are 
important when testing a full numerical implementation of these equations, important 
for generating the complete matter power spectrum for /(i?) gravity theories with a 
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realistic background cosmological expansion history. We also found that for i?" gravity 
to be consistent with the Meszaros effect, the parameter n needs to satisfy n > 2/3. 

We also gave a new covariant characterisation of the quasi-static approximation 
and used this to show that on small scales this approximation is valid for values of n 
in the neighbourhood of 1, i.e., it is in good agreement with a numerical integration 
of the full set of equations for the given set of initial conditions. This is the first time 
such a quasi-static analysis has been presented in a covariant way both for radiation 
and dust universes and provided the foundations for detailed comparison with what 
is found using the metric formalisms, together with a full computation of the power 
spectra. This will be presented in a future work. 
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